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This paper describes a simple numerical in- 
tegration method for diffraction integrals 
which is based on elementary geometrical 
considerations of the manner in which 
different portions of the incident wavefront 
contribute to the diffracted field. The 
method is applicable in a wide range of 
cases as the assumptions regarding the 
type of integral are minimal, and the results 
are accurate even when the wavefront is 
divided into only a relatively small number 
of summation elements. Higher accura- 
cies can be achieved by increasing the 
number of summation elements and/or 
incorporating Simpson's rule into the basic 
integration formula. The use of the 
method is illustrated by numerical exam- 



ples based on Fresnel's diffraction inte- 
grals for circular apertures and apertures 
bounded by infinite straight lines (slits, 
half planes). In the latter cases, the numeri- 
cal integration formula is reduced to a 
simple recursion formula, so that there is 
no need to perform repetitive summations 
for every point of the diffraction profile. 
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1. Basic Equations 

Diffraction problems that can be described in two 
dimensions as indicated in Fig. 1 usually lead to com- 
plex integrals of the type 



U{x,z) = S^^A{^,x,z)B{^,x,z), 



(1) 



where A and B account for the amplitude and phase of 
the optical field at a point of observation P = (jc,z) due 
to a point Q = (^,0) located inside a diffracting aperture. 
In practice, the phase term B oscillates rapidly inside the 
range of integration while the amplitude term A varies 
slowly. The phase term B is often, but not always, a 
sinusoidal function of the form e'''^'^^, where k = Iti/X is 
the circular wavenumber of monochromatic light with 
wavelength A . 



When closed analytical expressions for U are not 
available it is sometimes possible to find approximate 
solutions by the method of stationary phase [1,2]. How- 
ever, in this author's experience, the results obtained can 
be unreliable. A potential source of error is the basic 
premise of the method itself; namely, that the rapid 
oscillations of the phase term B nullify each other ex- 
cept in the vicinity of stationary points. This is pre- 
sumed true on account of the slow variation of the am- 
plitude term A, but inconsistent with the well-known 
fact that the corresponding series encountered in con- 
nection with Fresnel's zone construction, Ai — A2 
-I- A3 — ... — A„, has a finite value even when the terms 
alter their absolute values very gradually [3,4]. Accord- 
ingly, substantial errors can occur if the contribution of 
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Fig. 1. Cross section of a plane diffracting aperture (xz -plane). 



the Stationary points is weak. As shown in Appendix A, 
the stationary-phase method fails completely when 
applied to Fresnel diffraction at a slit or half plane. 

For these reasons, it may be preferable to use numer- 
ical integration. In doing so, the width of the summation 
elements, A^= /z, must be sufficiently small to ensure 
that all oscillations of the phase term B are accurately 
sampled. On the other hand, the computations would be 
unnecessarily complicated if h is too small. A general 
rule for choosing h can be established as follows. As we 
are dealing with diffraction, the period of the oscilla- 
tions does not depend on the specific functional form of 
B but only on the path length PoQ + QP, where 
Po = (xo,Zq) is the source point shown in Fig. 1. When the 
point Q is moved along the jc-axis by an increment h, 
this path length changes by 

S=Vz^ + (^+h- xof - Vzo' + (^ - -^o)' 



where it is assumed that (^ — Xq)^« zo^, H ~ x)^« z^, 
and terms in h^ are ignored. Hence it follows from the 
quarter-wave criterion that reliable results can be ex- 
pected when h is chosen so that 5 < A/4. 

We now divide the aperture halfwidth a into A^^ sum- 
mation elements bounded by equidistant points Q„, as 
illustrated in Figs. 2a, b, and 4, and define 



OQo = x = mh, QQQn=i-x, Qn-iQn = h = — , (3a) 



where O is the coordinate origin, Qo is the projection of 
the point of observation P onto the jc-axis, and m, n, and 
A'^ are integers. Accordingly, Eq. (1) can be replaced by 
the quadrature formula 






(3b) 



V-^ 



z' + i^+h-xf + ^z^ + i^ 



-xf 



where 



h\^-^^J 



Izol 



(2) 



A^n = A[{n — j)h, mh ,z], 5„„ = B[(n — ^)h, mh , z ] 

(3c) 
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(a) x<a 




(b) x>a 



Fig. 2. Annular summation elements for a circular aperture (xy -plane), (a) Lit region, (b) 
Shadow region. 



are the values of A and B at the midpoints of the summa- 
tion elements and the limits of summation depend on the 
diffraction problem being considered. The value of A'^ to 
be used in these formulae can be estimated from Eq. (2) 
by assuming a distant source (\zo\»z) so that 
8^h(^ - x)lz = a^nl{Nh). If this is to be less than A/4 
for every summation element used for the computations, 
a good upper limit for n is 3N ,^ and then one finds 



The value pertains to a point of observation located one aperture 
halfwidth beyond the shadow boundary. A larger value should be used 
if necessary. 



A^; 



12a- 
Xz 



6m 

— = 2m. 

TT 



(3d) 



where u = ka ^Iz is the familiar configuration parameter 
of Fresnel's diffraction theory for \zq\»z. As it is well 
known that diffraction patterns pertaining to large val- 
ues of M are highly structured [4], this result makes good 
sense in that it stipulates narrower summation elements 
when M increases. The accuracy of the numerical com- 
putations can also be improved by replacing the value of 
B,nn in Eq. (3c) with the value corresponding to Simp- 
son's rule. 



583 



Volume 105, Number 4, July-August 2000 

Journal of Research of the National Institute of Standards and Technology 



B^„ = ^[B[(n - l)h,mh,z] + 4B[(n -\)h,mh,z] 

+ B[nh,mh,z]}. (3e) 



given by 



dQ = 2iT(l-x/^)d(^-x)(^-x), 



(5b) 



By means of trial computations, it was found that this 
substitution can result in a tenfold improvement of accu- 
racy. 

The above equations are intended for applications 
where closed solutions of Eq. (1) cannot be found and 
will be used in future research. In the remainder of the 
present paper, their validity will be demonstrated by 
numerical examples involving the Fresnel diffraction 
integral 

f/F(P) = f/geom(P)«F(P), «f(P) = ^J^Q^'"^^'''^^ (4a) 

where t/geom is the geometrical field in the absence of 
diffraction, ctp is the modification of the field by diffrac- 
tion, and where 



QP = \/(i-xf + v' + z'-z + ^^ f^"^' . (4b) 

2.Z 



These expressions are valid in the paraxial Fresnel ap- 
proximation for a distant source point Po and pertain to 
a point of observation F = (x,0,z) as in Fig. 1 whereas 
the point Q = (^,17,0) is assumed to lie anywhere in the 
aperture plane z = 0. In Sees. 2 and 3 of the paper, ctp 
will be reduced to a two-dimensional integral for the 
respective cases of circular apertures and apertures 
bounded by infinite straight lines (slit and half plane). 
The results of the numerical integration will be shown 
and compared to the corresponding exact solutions, 
which may be found in Ref. [5]. 



2. Circular Aperture 

For a circular aperture of radius a, Fresnel's integral 
in Eq. (4a) can be reduced to a single integral by assum- 
ing annular area elements dQ which are centered on the 
projection Qo of the point of observation onto the aper- 
ture plane, as indicated in Figs. 2a and 2b. With O as the 
aperture center, x = OQo > and ^ = OQ > x for points 
to the right of Qo, the distance QP defined by Eq. (4b) 
will then be constant and equal to 



where x = ^AQoQ is the semi-angle subtended at Qo by 
the intersection of the area elements with the aperture 
rim, as indicated in the Figs. 2a and 2b. 

When Qo lies inside the aperture (x ^ a, as in Fig. 
2a), the innermost area elements with radii ^ — x 
< a — X are unobstructed and fully contained in the 
aperture (;t'= 0), and the outermost elements with radii 
^ — x> a-\- X are fully obstructed (;^ = tt). For the inter- 
mediate elements the angle x is found by applying the 
cosine theorem to the triangle OAQo shown in the fig- 
ure, so that 



cos;^': 



a^-x^-i^-xf 
2x{^-x) 



(5c) 



Thus, upon substitution of Eqs. (5a) and (5b) into Eq. 
(4a) and noting that 2Tr/(Az) = klz = u/a^. 



- m 
af{x) = — — 



&i$-x)i$-x)d 



ik(i-x)^l(2z) 



+ d(i - x)(i - x)(l - xMo 



iki^-xmiz) 



= [l-c 



iu(a-xy-l(2a^)-i 



IM 



d(i - x)(i - x)(l - xM^ 



ik{i-xfl{2z) 



(5d) 



where the first integral was reduced to an elementary 
expression by substituting \k{^ — x^llz as the integra- 
tion variable. When Qo lies outside the aperture {x> a, 
as in Fig. 2b), the inner elements with radii 
^ — X ^ X — a and the outer elements with radii 
^ — X ^ x + a are all fully obstructed (;^ = tt). In the 
intermediate region, Eq. (5c) applies once again^ and we 
have 



ap(x) = — i^ I di^ - x)i^ - x)il - xl^)^ 



ik(i-x)'l(2z) 



(5e) 



QP = z + {^-xfllz 



(5a) 



everywhere on dQ and the integration can be carried out 
over ^ alone. As the annular area elements are eccentric 
to the aperture they are, in general, partially obstructed 
by the aperture rim so that their effective areas will be 



The integrals in Eqs. (5d) and (5e) can now be identi- 
fied with Eq. (1), with 



Although X is discontinuous on crossing the shadow boundary, Eqs. 
(5d) and (5e) are equivalent and give identical results for x = a. 
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A{^,x,z) = 



luj^- x){l - xl'^) 



Table 1. The largest errors in diffraction profiles computed from 
Eqs. (6a), (6b), (6c), (10a), (10b), and (lib). 



A^^ 



accon 


B(^,x,z) 
ling to Eqs 


. (3a- 


c), 






(6a) 


N 


Circular 
aperture 
(m=5tt) 


Slit 
(m=5'tt) 


M 


Half 
plane 




16 


6.0X10"^ 


2.3X10"' 


8 


1.9X10"' 


«f(P) = 


= a„ = 


[1- 


iu(N-m)y(2N^)-i 






32 
64 


2.3X10"' 
6.4X10"' 


5.6X10"' 
1.4X10"' 


16 

32 


4.3X10"' 
1.1X10"' 


N+m 
















128 


2.5X10"' 


3.6X10"" 


64 


2.7X10"" 


n=N—m+ 


1 


~ Xmi 


y^)gi«(«-i)2/(2A^2)^ m 


<N 


(6b) 


256 
512 


9.6X10"" 
3.5X10"" 


9.0X10"' 

2.3X10"' 


128 

256 


6.6X10"' 
1.7X10"' 



«f(P) -«,„=- ^ E n(l- Xmn/^W""-^'""''"\ (6c) 

n=m—N 



where 



readily achieved by choosing cartesian coordinates so 
that the 3; -axis is parallel to the aperture edges. Letting 
F = (x,0,z) as before, this leads to 



Xfnn COS 



A^' - m' - (n - {f 
2m(n — 5) 



(6d) 



The use of these equations on a personal computer is 
simple. As an example, Fig. 3 compares the numerical 
values of la^P obtained for u = 5tt and A/^= 16 to the 
exact results given by the Fresnel-Lommel theory. The 
agreement is good and improves when larger values of 
N are used, as indicated in the left-hand column of Table 
1. The values listed in the table are the maximum errors 
encountered in the range m < 3A^. In this particular case 
they occurred near the center of the profile (m < O.SA'^). 




Circular aperture 



Fig. 3. Approximate ( — ) and exact ( — ) diffraction profiles of a 
circular aperture. 



3. Apertures Bounded by Infinite Straight 
Lines 

For a plane aperture of width (/ -i- r) which is 
bounded by infinite straight lines as in Fig. 4, the reduc- 
tion of the integral of Eq. (4a) to two dimensions is 



«P(P) = ^ I dTje'*^ '2z I ^^^m-.r,(2,) 



(r-x) 



le' 



VAz J 



(7a) 



where 



dTje 



ikT]H(2z) 



'2XzFi^) = VXzd 



(7b) 



and F(oo) = e''^''* is the complex Fresnel integral at infin- 
ity. 

The amplitude term A(^,jc,z) in Eq. (1) is now given 
by the factor — ie''^'V V Az that appears outside the inte- 
gral of Eq. (7b) so that, on letting 1 = Lh and r = Rh as 
indicated in Fig. 4 and using Eqs. (3a) to (3d), we find 



ihe' 



Az „=- 



E 



,ik{n-^)^h^/(2z) 



n^O. (8a) 



(L+m) 



It follows immediately that, if a„, is known and P is 
moved to the right or left hy ±h, the new value of a,„ 
is given by the recursion formula 



^m±l ^m ~l~ 



i"^ r^ik(L+m+j)^h^/(2z) _ ^ik(R-m-^)hM2zh /g^\ 



Az 



which illustrates in a very instructive manner how the 
diffraction pattern changes when the point of observa- 
tion is moved, so that one portion of the wavefront is 
uncovered by the aperture and another portion is cov- 
ered. The recursion formula Eq. (8b) is convenient for 
practical computations as it allows the computation of 
successive values of a,„ without performing the summa- 
tion of Eq. (8a) for every point. In the examples given 
below, a starting value for a„ is obtained from known 
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Fig. 4. Rectangular summation elements for apertures bounded by straight lines (xy -plane). 



closed solutions, so that there is no need at all to perform 
a summation. This use of an exact starting value also 
improves the accuracy of the computations because it 
forces the initial error to be zero. 



3.1 Slit 

For a slit of width 2a , we define 1= r= a, L = R = N 
and h = a/N, so that 



h 



1 



kh' 



V\z A^V27r' z 



u 



(9a) 



X -- + - 

'" ~ 4 2 



1 + 



(m + {) 



N' 



Y = 



u(m + ^) 

N 



(10c) 



As an example, Fig. 5 shows the approximate and exact 
diffraction profiles of a slit for m = 5tt, A'^= 16, and 
using the well-known Fresnel solution, aF(0) 
= (1 — /)F( V m/tt), as the starting value. The two 
curves resemble each other closely, the largest errors 
being on the order of 0.023 and occurring near 
mIN = 0.5. The improvement of accuracy achieved by 
using larger values of A'^ is shown in the center portion 
of Table 1. 



^ _ ^ ^,J:L. plTTMr iH(A'+m+i)2/(2/V2) _ iH(A'-m-l)2/(2/V2)-| 



(9b) 



where it is noted that it is sufficient to perform the 
computations for m > as the diffraction pattern is sym- 
metrical. 

On account of the trigonometric identity 



.i/3_^i(a+/3)/2rgi(a-/3)/2 _ ^-i{a-p)l2^ 



e''^ = e" 



^2ie'(«+/3V2sin[(a-/3)/2], 



(10a) 



this can be transformed into the following expressions, 
which are convenient for practical computations: 



_ J_ /2m -^x^ ■ y 

N ^ TT 



(10b) 




Slit 



Fig. 5. Approximate ( — ) and exact ( — ) diffraction profiles of a slit. 
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3.2 Half Plane 

On letting L = and i? = oo, Eqs. (8a) and (8b) apply 
to a diffracting straight edge which coincides with the 
3; -axis of Fig. 4. As there is no aperture edge on the 
right, the last term of Eq. (8b) is now absent and one 
obtains 



where (^=^-;c,A = - ie^'VVAz, gicf)) = (j)'lz. The 
only stationary point, defined by g'(4>) = 4>lz = 0, oc- 
curs at (/)o = so that, according to the Appendix of Ref . 
[2], the stationary-phase value of the integral is 



a^ix) ~ 



Va^ 



,i«r(m±l-i)2/i2/(2z) 



Az 



(11a) 



The previous definition of /z as a given fraction of 
aperture width is no longer applicable but can be 
replaced by taking /z as a certain fraction, say h 
= vXz/M, of the width of the first Fresnel zone. Hence 
it follows easily that 



"" "^M ^ 



X„ 



1 , (m ± 1 - ^f 



V 



M' 



(lib) 



Figure 6 shows the diffraction pattern computed in 
this manner for M = 4, using the well-known Fresnel 
solution, aF(0) = V2, as the starting value. The agree- 
ment with the exact result is within ± 0.025 for 
— 4 < w/M < 1, but considerably worse beyond these 
limits. These discrepancies are reduced when larger val- 
ues of M are used, as may be seen from the right-hand 
side of Table 1 . 




1.5 m/M 



Fig. 6. Approximate ( — ) and exact ( — ) diffraction profiles of a half 
plane. 



4. Appendix A. Evaluation of Eq. (7a) by 
the Method of Stationary Phase 



Vl^"(cAo)l 



j^ QikgCpo) 



where ^"((^o) = 1/z and s = e"^™'"* according as g"((l)o) is 
positive or negative. Since 1/z is positive, this gives 
ap~ I so that the "diffraction pattern" of a slit or half 
plane would be identically equal to one everywhere in 
the plane of observation. 
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Equation (7a) is of the form 



(r-x) 

aj,ix)= I dc^Ae'*^^'^^ 



-il-x) 
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